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ABSTRACT 

It  has  been  previously  shown  that  the  temporal  integration  of  hyperbolic  partial  differ¬ 
ential  equations  may  ,  because  of  boundary  conditions,  lead  to  deterioration  of  accuracy  of 
the  solution.  A  procedure  for  removal  of  this  error  in  the  linear  case  has  been  established 
previously. 

In  the  present  paper  we  consider  hyperbolic  p.d.e’s  (linear  and  non-linear)  whose  bound¬ 
ary  treatment  is  done  via  the  SAT-procedure.  A  methodology  is  present  for  recovery  of  the 
full  order  of  accuracy,  and  has  been  applied  to  the  case  of  a  4th  order  explicit  finite  difference 
scheme. 
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1.  Introduction 

Currently  there  is  a  growing  interest  in  long  time  integration  for  solving  problems  in  areas 
such  a  aero-acoustics,  electro-magnetics,  material-science  and  others.  This  necessitates,  ([!]), 
working  with  higher  order  (4th  order  accuracy  and  above)  schemes.  Often  the  methodology 
of  choice  is  to  semi-discretize  the  equations  by  applying  a  high  order  (4th  and  above)  spatial 
difference  operator  and  then  advance  temporally  using  single  level  multi-stage  Runge-Kutta 
integrators.  This  raises  the  question  of  how  to  supply  boundary  values  at  the  intermediate 
stages  of  the  R-K  integrations.  In  the  case  of  hyperbolic  p.d.e’s  that  means  the  imposition 
of  the  time  dependent  conditions  at  the  inflow  boundary. 

The  conventional  (and  intuitively  natural)  way  of  imposing  inflow  boundary  conditions 
at  the  intermediate  stages  is  to  use  the  “appropriate"  value  of  the  boundary  data,  g{t)^  at 
each  stage.  Thus,  for  example,  at  a  stage  corresponding  to  t  -f-  (Ai/2),  one  would  impose 
g{t  -f  A</2). 

In  a  previous  paper  ([2])  ,  it  was  shown  that  the  procedure  described  above,  when  ap¬ 
plied  to  hyperbolic  p.d.e’s  with  time  dependent  b.c.’s,  reduces  the  accuracy  near  the  inflow 
boundary  to  first  order  and  thus  the  overall  accuracy  cannot  exceed  0(A.x^).  This  conclusion 
is  independent  of  the  order  of  accuracy  of  the  spatial  difference  operator. 

One  way  of  avoiding  the  dilemma  of  what  boundary  values  one  should  supply  at  the 
intermediate  stages,  is  to  advance  the  R-K  integration  without  imposing  any  intermediate 
values,  but  rather  obtain  the  intermediate  boundary  values  from  the  numerical  solution 
operator.  However,  this  approach  has  the  disadvantage  of  reducing  substantially  the  stability 
limit  (e.g.,  the  allowable  time  step  is  reduced  by  a  factor  2  in  the  case  of  4th  order  classical 
R.-K  with  a  4th  order  spatial  derivative  operator),  hence  rendering  it  less  than  attractive. 

In  ([2])  a  general  methodology  was  presented  in  the  case  of  linear  p.d.e.’s  for  the  correct 
imposition  of  the  intermediate  stage  boundary  values  so  that  the  scheme  recovers  its  full 
formal  accuracy.  This  was  expounded  in  detail  for  the  case  of  the  classic  4th  order  R-K 
with  a  hyperbolic  4th  order  spatial  difference  operator.  It  was  also  shown  there  that  in  the 
non-linear  case  (e.g.  hyperbolic  conservation  laws)  this  methodology  was  applicable  to  R-K 
integration  up  to  3rd  order  .  For  R-K  .  methods  of  4th  order  and  above  we  were  not  able  to 
extend  the  theoretical  approach  described  in  (  [2]). 

In  this  paper  we  address  anew  the  issue  of  how  to  deal  with  the  non-linear  case.  We 
present  a  methodology  for  retaining  the  full  accuracy  even  in  the  non-linear  case.  The 
application  of  this  methodology  involved  numerical  determination  of  free  parameters  in  con- 
tradistinctsr  numerical  determination  of  free  parameters  in  contradistinction  to  the  linear 
procedure  described  in  ([2]).  We  And  for  example,  that  in  the  4th  order  classical  R-K  integra¬ 
tor  with  4th  order  explicit  spatial  derivative  operator,  the  full  accuracy  is  retained  without 
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any  reduction  in  t  he  allowable  time  step. 

The  new  procedure  is  demonstrated  for  hyperbolic  problems  where  the  boundary  condi¬ 
tions  are  satisfied  by  using  the  SAT  approach  ([3].  The  motivation  for  doing  so  is  that  the 
SAT  procedure  is  the  only  one  that  prevents  temporal  growth  not  present  in  the  true  solu¬ 
tion  of  a  system  of  p.d.e.’s.  Section  2  describes  how  to  apply  correctly  the  intermediate  SAT 
boundary  conditions  in  the  case  of  a  linear  problem.  In  Section  3  we  cover  the  non-linear 
case. 

2.  The  Linear  Case 

In  this  section  we  analyze  the  effect  of  imposing  the  inflow  boundary  conditions  in  the 
conventional  way  when  the  discretization  algorithm  employs  the  SAT  approach  (see  ([3])  ). 
The  SAT  is  a  penalty  type  method  that  was  constructed  so  as  to  ensure  that  the  numerical 
solution  will  not  include  temporal  growth  which  is  not  of  a  physical  origin.  This  is  achieved 
by  mimicking  the  energy  estimate  of  the  p.d.e. 

Recall  that  we  are  considering  the  following  hyperbolic  problem  (see  ([3])): 

=  0<.T</,  <>0  (2.1) 

u(0,t)=ir(0  (2.2). 

The  SAT  formulation  for  the  semi-discrete  version  of  (2.1)  -  (2.2),  based  on  a  uniform  grid, 
is 

=  {jO  V  it))  -  (vo(i)  -  g{i)),  !  =  (>0  (2.3) 

where  V=  [fo, •  •  • , ^’yv]^  is  the  semi-discrete  approximation  that  converges  to  u(x,-, f)  at 
the  spatial  grid  points  .t,-  (for  stable  discretizations);  and  j-^D  is  the  differential  matrix  rep¬ 
resentation  of  the  derivative  operator  (— ^)-  The  vector  a  depends  on  the  differentiation 
matrix  ^D,  and  on  the  energy  norm  used  in  bounding  the  error.  It  is  determined  as  described 
in  ([3]);  see  the  discussion  after  equation  (6)  therein. 

The  demonstration  of  accuracy  deterioration  will  be  shown  for  the  four  stage  “classical” 
RK  algorithm,  which  is  one  of  the  most  commonly  used  RK  time  advancing  schemes.  For 
the  analysis  to  make  sense  we  assume  that  the  spatial  discretization  is  at  least  fourth  order 
accurate. 

The  above  mentioned  four  stage  RK  integrator  is  implemented  as  follows: 

1/(1)  =  ]/(")  +  ^  [4")  _  g(t)]  (2.4) 
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(2.5) 


y(2)  =  y(n)  J^(i)  _  ^ 

y(3)  ^  y{n)  ^  XDV^^^  -Xra  4^^  "  ^  0  "'■  y) 

y(n+i)  _  y{n)  ^  A X)[y(")  +  2V(^)  + 
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-^T^{[^'o"^-<7(0]  +  2[4‘^-<if(^  + y)  +2  + 

+[v^q'‘ -  git  +  ^i)]}  ■ 

To  check  for  accuracy  we  substitute  for  the  the  exact  values  ii{xi,t),  and  in  particular 
=  g{t).  Note  from  eq.  (2.4),  that  on  the  boundary  (using  the  differential  equations 
(2.1), (2.2)) 

4^)  =  4")  _  f  =  g[t)  +  ^g\t).  (2.8) 

Thus  in  eq.  (2.7),  we  have  for  the  term  Vq^  -  g  ^i  + 

4'’  -  s  ('  +  t)  =  [«(')  +  “  «  {'  +  f)  = 

Thus  -  yf")  is  at  best  0(A<^),  and  not  0(At®)  as  expected  from  the  R-K  scheme  used. 

In  this  linear  case,  the  remedy  proposed  in  the  previous  paper,  ([2]),  works  here  as  well. 
In  particular,  eqs.  (2.4)  -(2.7)  take  the  following  form: 


1/(1)  =  1/("1  +  jDl/(">  -^TC  -  s(()l 

(2.10) 

1/(2)  =  l/(»)  +  Aci/d)  _  ?  [«*'>  -  g{t)  -  y5'(() 

(2.11) 

Ai  A/2 

y(3)  ^  y(n)  ^  ^£,^(2)  _Xra  4"'^  -  g{t)  -  ^4(0  “  -J-9"i^) 

(2.12) 

yn+l  ^  yn  Jy(n)  _|.  -f 

-  Y<t  |[4”^  -  git)]  -1-  2  4*^  -  git)  - 

+2 1^4  ^  -  git)  -  2  ^  4  ^ 

(2.13) 

+  4  -  g(t)  -  ^tg  - 

It  is  readily  verified  that  1/1"+^)  —  f/f"!  =  0(A<®),  as  required. 
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3.  The  Non-Linear  Case 

For  the  sake  of  simplicity  we  consider  first  the  scalar  conservation  law  p.d.e. 

^  +  =  0<a;<l;  t>0  (3.1) 

u{0,t)=g{t).  (3.2) 

In  general,  for  any  spatial  discretization  (whether  explicit  or  implicit)  the  semi-discrete 
form  of  (3.1)  -  (3.2)  is: 

^  =  jDUV)  -  ir  ?  -  9(()1.  (3.3) 

Using  the  notation  of  reference  ([3]), 

D  =  -P-^Q  (3.4) 

where  is  the  differentiation  matrix  representing  the  differential  operator,  (-^);  com¬ 
posed  of  the  explicit  part  Q  and  the  inverse  of  the  implicit  part  P.  For  a  fully  explicit  spatial 
differentiation,  P  =  I  +  B,  where  B  differs  from  zero  only  at  the  two  diagonal  corners  (see 
examples  of  P  and  in  ([4])).. 

The  vector  cr,  again  using  the  notation  of  [  ],  is  given  by 

a=  ha{uo)gooP~^H-'^S  (3.5) 

where  o(uo)  =  (|f)  and  5^00  is  twice  the  value  of  the  left  upper  corner  element  of  HQ. 
For  the  definition  of  the  matrix  H  see  Assumption  I  in  [1].  The  parameter  r  is  determined 
from  stability  consideration  to  be  r  >  1,  see  ([3]). 

Next  we  demonstrate  that  writing  the  classical  4th  order  RK  for  eq.  (3.3),  using  the 
linear  “fix”  as  in  equations  (2.10)  -  (2.13),  does  not  yield  the  required  4th  order  accuracy: 

,/(»  =  y(-)  +  -  jT  a  (4”'  -  9(01  (3.5) 

VW  =  yl")  +  -  jT  5  -  g{t)  -  y9'(0]  (3.6) 

A/  AP 

y(3)  =  ,/(.)  +  XDfiV^^')  -\ra  -  9(0  -  “j'lO  -  — 9"(0  (3.T) 
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,/(»+!)  =  v'l"'  +  j-Dl/(V")  +  2f(V''>)  +  2f(V'^')  +  /(V<=')| 
-^TT  -  g(t)]  +  2  [4"  -  g(t)  -  ys'(()] 


+2U^^-9(t)-fg'(t)-^g"{t) 


4  -  g{t)  -  -  -y9  -  ~9  I  • 


Again,  when  checking  accuracy,  we  take  =  «(rc,-,  t),  and  in  particular  Uq  ^  —  9{i)] 
also  \DJ{V'^)  =  -^/(w)  +  0(At‘').  With  these  preliminaries  we  get  immediately  from 
equation  (3.5) 

4"  = 

Note  that  this  is  the  same  as  in  the  linear  case,  see  eq.  (2.8).  Thus  with  Vq  ^  =  g{t)  and 
with  equation  (3.9)  we  can  supply,  for  the  purpose  of  accuracy  checking,  the  correct  values 
of  and  When  we  look  at  eq.  (3.6),  using  the  above  results,  the  governing  non-linear 
p.d.e.,  and  simple  Taylor’s  expansion,  we  have: 

i/(2)  =  +  y(n)  + 

=  -  ytJ  ['■  -  f  h  f  ¥ + 

=  u  — —  —ut  — —utt  +  O(At^). 

(3.10) 

So  finally  we  have  on  the  boundary 

t,W=5(()  +  Ys'(l)  +  ^/(()  +  0(A(=’),  (3.11) 

It  follows  from  (3.11)  that  the  “penalty”  term  in  (3.8)  introduces  an  error  of  (At^).  Since 
the  coefficient  of  (A/.^),  [/„„(w<)^]o»  cannot  generally  be  expressed  as  a  function  of  g{t)  and 


(3.11) 
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its  derivatives,  it  is  very  difficult  to  remedy  the  situation.  Thus  the  “linear  procedure”  fails 
at  the  third  RK  stage. 

We  now  propose  a  methodology  to  deal  with  the  RK  integration  of  non-linear  hyperbolic 
conservation  laws.  We  first  present  this  procedure  in  the  case  of  the  classic  4th  order  R-K 
scheme.  Our  starting  point  is  the  observation  that  the  “linear-procedure”  yields  the  required 
accuracy  for  ‘‘^nd  The  idea  is  to  use  at  each  stage  a  linear  combination  of  the  “linear” 
SAT,  or  penalty,  terms  used  (3.5)  and  (3.6).  The  4th  order  classic  RK  stages  will  thus  be; 

yd)  =  -f-  ^Z)/(K”)  (3.12) 


7(2)  =  7(")  +  ^J[l/(l/(i))  [4^^  -  g(t)  -  ^g'it)]  -  ^7  ^  [4"^  -  9{t)]  (3-13) 

7(3)  ^  y(n)  ^  XDfiV^^^)  -X8a  [4")  -  g{t)]  -Xea  [vi^^  -  g{t)  -  ^g'{t)]  (3.14) 

yin+i)  ^  yin)  ^  ^D[/(7")  +  2/(l/(*))  +  2 /(l/(2))  -h  /(!/("))] 

(3.15) 

-pt  ^  bo”’  -  <i(0]  -  ^  4’’  -  9{t)  - 

where  the  free  parameters  a,  /3, 7, 6,  e,  /i,  and  1/  will  be  chosen  so  as  to  maximize  the  allowable 
time  step.  It  is  clear  from  the  previous  discussion  that  the  system  (3.12)  -  (3.15)  maintains 
the  4th  order  accuracy.  There  remains  the  question  of  whether  the  CFL  stability  condition 
deteriorates,  in  comparison  to  the  conventional  application  of  R-K,  equation  (2.4)  -  (2.7).  It 
is  also  clear  that  the  optimal  choice  of  the  free  parameters  a, . . . ,  varies  with  the  spatial 
discretizations  (i.e.  the  differentiation  matrix  D)  and  boundary  closures.  One  checks  that 
the  absolute  values  of  the  eigenvalues  of  the  amplification  matrix  resulting  from  (3.12)  - 
(3.15)  should  not  exceed  unity.  We  carried  out  this  procedure  (using  Matematica  Software) 
in  the  case  of  an  explicit  4th  order  algorithm  with  3rd  order  boundary  closures.  In  this  case 
H  =  /,  and  the  matrices  P  and  Q  can  be  found  in  section  (9.1)  of  reference  [3].  With  the 
following  values  of  the  free  parameters,  a  =  ^  =  — £  =  1,  <5  =  2,  /<  =  0,  =  3,  7  =  —0.37, 

the  CFL  condition  becomes  A  <  2.1.  This  is  the  same  restriction  on  the  time  step,  A^,  as 
one  has  in  the  linear  problem  using  the  “conventional”,  i.e.  less  accurate,  boundary  values 
with  or  without  the  SAT-term. 

4.  Conclusions 


In  summary,  we  have  in  (3.12)-(3.15)  a  4th  order  RK  scheme,  applied  to  a  non-linear 
p.d.e.,  which  maintains  the  overall  4th  order  accuracy  without  any  decrease  in  the  allowable 
time  step.  The  extension  to  a  system  of  hyperbolic  p.d.e. ’s  is  quite  straight  forward  using 
the  SAT-system  approach  delineated  in  ([3]). 
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